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1  INTRODUCTION 


As  with  Fourier  analysis,  the  basic  idea  behind  wavelet  analysis  is  to  decompose  a  function  into  the 
sum  of  meaningful  'basis'  functions.  Unlike  Fourier  methods,  however,  the  basis  functions  of 
wavelet  analysis  exhibit  some  form  of  localisation  in  both  the  frequency  and  time  domains.  The 
desirability  of  such  a  decomposition  is  well  recognised  and  has  led  to  the  development  of  the  short 
time  Fourier  transform.  In  what  follows,  it  will  be  shown  that  wavelet  techniques  have  several 
features  that  make  them  an  attractive  alternative  to  the  short  time  transform. 

The  following  report  describes  some  of  the  basic  ideas  that  underlie  a  discrete  wavelet  analysis.  A 
non  rigorous  approach  is  used,  but  there  are  several  references  to  more  complete  descriptions.  The 
major  purpose  is  to  give  the  reader  sufficient  understanding  to  allow  a  meaningful  application  of  the 
software  that  was  developed  as  part  of  this  project  Appendix  A  contains  some  FORTRAN 
subroutines  that  are  sufficient  for  a  basic  wavelet  analysis.  In  addition.  Appendix  B  describes  an 
associated  PC  software  package  that  can  perform  a  fairly  comprehensive  analysis.  Several  examples 
of  its  use  are  included  and  these  amply  demonstrate  the  potential  of  a  wavelet  approach. 

The  report  draws  heavily  on  some  papers  by  Daubechies  [1]  and  Mallat  [2,3].  In  addition,  extensive 
use  has  been  made  of  an  excellent  introductory  paper  by  Strang  [4]  and  a  recent  tutorial  paper  by 
Rioul  and  Vetterli  [5].  For  more  advanced  material,  the  reader  should  consult  the  books  by 
Daubechies  [6]  and  Ruskai  et  al  [7]. 

2  THE  SCALING  FUNCTION 


In  Fourier  analysis,  the  basis  functions  are  generated  by  scaling  the  argument  of  a  single  function  ({) 
(({)(0=exp(fr)).  If  the  sequence  <t)^  («eZ)  is  generated  according  to  (1)^(/)  =())(«/),  it  forms  a  basis 

for  (the  space  of  Lebesgue  square  integrable  functions  on  [-7t,7t]).  In  the  case  of 

wavelet  analysis,  the  basis  elements  are  generated  from  a  single  function  by  both  a  scaling  and  a 
translation.  Consider  the  problem  of  representing  a  function  in  terms  of  elements  (|)|^  iiJeZ)  that 

j 

are  generated  from  the  single  function  <()  according  to  ({)y=2^4»(2^x-0.  The  simplest  example  of(t) 
is  the  box  function 


<t>W  = 


1  0<x<l 

0  otherwise 


which  leads  to  a  step  approximation  A  J  for  the  function  /  where 


(1) 


(2) 


and  the  value  of  J  sets  the  level  of  approximation.  As  illustrated  in  Figure  1,  the  accuracy  will 
inaease  as  J  grows  {2'^  is  the  width  of  a  box  at  level  J  if  the  boxes  have  unit  width  at  level 
7=0). 

For  the  above  ((> ,  it  is  of  interest  to  note  that 

((>(jt)  =  <t>(2jt)  +  <))(2x-l) 
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In  general,  il  will  be  required  that  satisfies  the  dilation  equation 

<t»(jc)  =  53  c\i^2x-k) 

k 


(4) 


a  property  that  is  crucial  to  what  follows.  Up  to  a  suitable  normalisation  factor,  relation  (4)  will 
define  (().  The  iiormalisation  =  1  is  chosen  for  the  current  considerations  and  this  imposes 

the  condition  53  ^  t  ~  ^  •  Function  <()  is  known  as  the  scaling  function. 


If  coefficients  are  zero,  except  possibly  for  some  k  between  0  and  N,  it  can  be  shown  [1]  that 
the  scaling  function  will  be  zero  for  points  outside  the  interval  [0,N].  (Except  for  the  box  function, 
it  will  nonnally  be  required  that  <|)(0)  =(t>(AO=0.) 

From  (4),  the  Fourier  transform  of  (t>  is  given  by 

=  ^CjJ<!)(2ac-l:)exp(r(0ix)dx  (5) 


and,  on  setting  y~2x-k. 


^((0)  =  J<t)(y)expfi^\fy  I  '53cjexp[ 


y 


•  ■  I  2 


Define 


2  * 


then 


(^((0)  =  (J) 


J. 


0) 


from  which 


^(to)  =  n  ^ 


^20 


(0 


2n 

\  ) 


In  the  limit  N—^,  this  implies 


4  =  11'’ 

j-i 


f  0)^ 

yj 

j 


(6) 


(7) 


(8) 


(9) 


(10) 


since  ^(0)  =  1 . 
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Condition  A  P  has  a  zero  of  order  p  ai  (o  =  rt . 

For  m=0,l,...,p-l,  this  implies  that 

53(-l)‘*"c-^  =  0  (11) 

k 

If  condition  A  holds  and  f  is  a  smooth  function,  coefficients  can  be  found  [6]  such  that 

I/-E  ^  ( 12) 

k 

for  some  constant  C.  Such  a  result  justifies  the  use  of  (2)  as  an  approximation  to  / . 


3  CALCULATION  OF  THE  SCALING  FUNCTION 


Consider  relation  (4)  at  integer  points  from  1  to  N-1  (c^=0  for  k>N  and  k<0)  .  This  provides 
the  system 


where  7"  has  a  maximum  eigenvalue  of  1.  The  power  method  can  be  used  to  find  the  corresponding 
eigenvector  and,  up  to  a  scaling  factor,  the  components  of  this  eigenvector  will  provide  values  for 
(!)(l),...,<t)((V-l).  Relation  (4)  can  be  used  to  find  intermediate  values  and,  after  sufficient  of  these 

have  been  calculated,  the  normalisation  imposed  by  use  of  a  suitable  quadrature 

procedure. 

Condition  O 

or 

(15) 

k 

For  fixed  J,  this  implies  that  the  <t)_y  are  ortbonormal. 
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Conditions  O  and  A  are  desirable  properties  and  will  be  satisfied  [  1  ]  if  P  has  the  form 

nm)  -  f (1‘) 


where 


L-\ 


\Q(a»\'^  =  E  j 


L-l*j 


,  .  1  (0.,  /  •  1  to,,  COSO), 

(sin- y  +  (sin* — //?( _ ) 

2  ~  ~ 


(17) 


and  R  is  a  real  odd  function. 

The  box  scaling  function  corresponds  to  L  =  I  and  ^=0  so  that  P((i))  =-l — Although  this 

scaling  function  is  discontinuous,  others  of  the  same  class  become  smoother  as  L  increases  (at  least  C 
for  L-l  ,  C'  for  L=4  and  C"  for  L=7). 


4  THE  WAVELET  FUNCTION 


Consider  the  approximations  Aj^J  and  A jf.  It  can  be  shown  [6]  that 

A,J  -AJ^DJ 


where 


(18) 


(19) 


and  V|;_^  =2^\li(2^;c-0  for  a  suitable  wavelet  function  \jr.  In  the  case  of  the  box  scale  function,  the 
corresponding  t|/  is  given  by 


1  ice  (0,^) 

^  -1  ice(_,l) 

2 

0  otherwise 


and  is  known  as  the  Haar  wavelet.  It  will  be  noted  that 

\|r(j:)  =  <t)(2x)-<{)(2jt-l) 

In  general,  the  wavelet  corresponding  to  the  scale  function  ((>  is  given  by 

¥(-«)  =  J^(-l)*c,.4(t>(2x-^) 


(20) 


(21) 


(22) 
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This  construction  ensures  that  t|/(x)  orthogonal  to  $(x-m)  for  meZ.  Furthermore,  comliuon  O 
implies  that  (i&Z)  forms  an  orthonormal  basis  for  the  Dj'.  Repeaung  the  reducuon  of 
equation  (18) 

AjJ  =  (23) 


and  taking  limits  in  J  and  K, 

/  =  E  E//  V, 


(24) 


It  can  be  shown  [  1 )  that  ( i,j&  Z)  form  a  orthonormal  basis  for  functions  in  L  '(R)  (the  space  of 
Lebesgue  square  integrable  functions). 

If  condition  A  holds,  wavelet  functions  have  the  following  important  properties: 

(1)  jx  ”'\^(x)dx=0  for  m=0  to  p-1. 

(2)  For  a  C”  function  /,  d,  <C2  -  . 


Daubechies  [1]  has  studied  scaling  functions  that  are  defmed  by  (16)  and  has  calculated  the 
coefficients  for  R-^  and  L  from  1  to  10.  For  example,  the  box  function  corresponds  toL=l 

and  has  Cq=c,  =  1.  For  L=2,  c,=(3V^V4,  and 

c^=0  otherwise.  Figures  2  ,3  and  4  show  the  scaling  functions,  wavelet  functions  and  their 
corresponding  Fourier  transforms  for  the  cases  where  L=l,  2  and  7.  From  these  examples,  it  can  be 
seen  that  the  scaling  function  has  the  nature  of  a  low  pass  filter  and  the  wavelet  function  the  nature 
of  a  band  pass  filter. 


5  WAVELET  ANALYSIS 


For  most  applications,  the  function  of  interest  is  only  known  in  terms  of  a  finite  number  of  samples. 
In  order  to  perform  a  wavelet  analysis,  these  samples  must  be  turned  into  an  initial  approximation 

Ajf.  A  wavelet  analysis  decomposes  Ajf  into  its  constituent  projections 

Dj  ^f,  Dj  J,  ... ).  Equivalently,  the  coefficients  {a/}  of  Ajf  are  decomposed  into  the 

coefficients  ({d/  'j,  [d-'^),  ...)  of  the  constituent  projections.  A  major  advantage  of  wavelet 
analysis  is  the  existence  of  fast  algorithms  for  performing  this  reduction.  The  present  section 
describes  these  algorithms  and  their  application.  In  addition,  the  interpretation  of  a  wavelet  analysis 
is  discussed. 
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Consider  the  decomposition 

or,  in  terms  of  the  basis  functions. 


(25) 


(26) 


where  a/  and  d/  .  (The  inner  product  <f,g>  is  defined  by  j/{u)g(u)du.) 

R 

Given  what  are  a/  and  d/ From  the  definition  of  (t>_^  and  relation  (4), 

a/  =  </■<!>, y> 

J 

-  2'^  jj[u)(f {2^ u  -i)du 

R 

=  2''^Cjj}(M)<{>(2-"‘M-2i-<:)dM 

=  2^52  ^K-.. 

R 


(27) 


if  K=2t +<:.  Since  ai^''  jf,u)^(2'''u-K)du,  it  follows  that 


/  1 

a.  =  -F=E  ^k-2.  « 

y/2  < 


/♦I 


(28) 


Consequently,  to  form  the  coefficients  [a/],  convolve  coefficients  (a/*'}  with  the  filter 
c  ^ 

coefficients  { - }  and  keep  alternate  samples  (a  low  pass  filtering).  In  a  similar  fashion 

<12 


■i.'  -  “ 

yfl  ^ 


y*i 

K 


(29) 


Consequently,  to  form  the  coefficients  [d/],  convolve  coefficients  {a/*’}  with  the  filter 

(-l)V  ... 

coefficients  { _ 1  and  keep  alternate  samples  (a  band  pass  filtering). 

(/2 
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Ii  IS  also  possible  to  reconstruct  the  civtficienis  |ti/ ')  trom  the  coefficients  {tj/l  and  [d/] 
From  (2X),  it  follow^  that 


(30» 


and,  from  (29), 


'  '  I  -« *Zt 


(31) 


Consequently,  on  taking  the  iniier  product  of  with  (26), 


j-i 


—Vu/c ,  -  _L5:tf/(-i)‘c, 
/1 1  ^ 


(32) 


In  effect,  coefficients  { «/*’ )  are  formed  from  coefficients  {a/  \  and  { d/ )  by  adding  the  results  of 
two  convolutions.  The  first  convolution  involves  the  coefficients  {«,  )  and  (  — and  the  secoid 


involves  the  coefficients  {d/ }  and  { 


In  order  to  produce  a  practical  implementation  of  the  above  procedures,  it  is  necessary  to  work  with 
a  finite  amount  of  data.  Unfortunately,  the  decomposibon  of  ff'*'  =(ao'^*',a/*’,...,tZ;5*J,)  into 

of  =  {afj\a/,...,asii )  and  )  requires  elements  that  are  not  contained  in  o'*’ .  If 

these  elements  are  unknown,  one  possible  procedure  is  to  zero  pad  the  known  elements.  It  should  be 
noted,  however,  that  this  approach  will  lead  to  an  increasing  number  of  incorreci  coefficients  as  the 
decomposibon  proceeds.  After  several  steps  of  decomposibon,  many  of  the  highest  N,  and  lowest 

N,  components  of  the  £■'  and  d^  will  be  contaminated.  An  altemabve  approach  is  to  assume  that 
the  inibal  data  represents  one  period  of  a  periodic  funebon.  It  should  be  noted,  however,  that  a 
periodic  funebon  is  not  square  integrable  and  so  a  complete  decomposibon  into  wavelets  is 
impossible.  In  this  case,  the  decomposibon  will  be  meaningless  beyond  K  steps  for  inibal  data 
consisbng  of  2^  samples. 
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A  further  problem  in  the  analysis  of  discrete  data  is  the  derivation  of  the  initial  coefficients  { a/ } . 
These  ;ire  calculated  according  to 


I 


=  <fAi> 

=  ^u)^(l'ii-i)du 

R 


(33) 


In  effect,  the  coefficients  are  obtained  by  sampling  the  function  after  passing  it  through  a  low  pass 
filter  that  corresponds  to  a  suitably  scaled  <().  Consider  the  samples  {//}  taken  from  /  at  intervals 
of  2’^.  The  coefficients  (ti/}  can  be  approximated  according  to 


J  i*/V 

«/  -  2'^  Y.fa-^kmk-i) 

km 


or,  in  a  more  convenient  form. 

-  2 


(34) 


(35) 


where  <t)^=(j)(fc).  This  is  an  approximate  implementation  of  the  filtration  and  sampling  processes. 
Appendix  A  contains  some  FORTRAN  subroutines  for  implementing  (35)  in  the  case  of  periodic 
data  that  has  been  sampled  at  unit  intervals  (initial  J  zero). 

The  wavelet  decomposition  is  economical  in  terms  of  both  storage  and  computation?)  requiiemenL 
Firstly,  the  decomposition  algorithm  only  requires  0{M)  operations  if  the  origina  data  set  hasM 
elements.  Secondly,  the  data  storage  can  be  accomplished  in  the  same  area  as  the  original  data  sincea'' 

becomes  which  then  becomes  and  so  on.  Appendix  A  contains  some 

FORTRAN  subroutines  for  the  wavelet  decomposition  (and  reconstructi'^n)  of  periodic  data. 

A  wavelet  decomposition,  can  be  crudely  regarded  as  a  spectral  analysis  with  frequency  space 
divided  into  bands  of  width  proportional  to  their  centre  frequency  (constant  Q  filtration). 

Furthermore,  the  analysis  exhibits  some  localisation  in  time  (the  domain  of  /).  Referring  to 
Figure  5,  the  basis  function  corresponding  to  a  given  wavelet  coefficient  will  make  its  major 
conuibution  over  the  region  that  is  delimited  by  the  box  surrounding  the  coefiicient  Crudely,  the 
square  moduli  of  the  wavelet  coefficients  yield  the  distribution  of  energy  in  terms  of  time  and 
frequency.  Referring  to  Figure  6,  the  various  projections  {Djf,  Dj_J,  ...)  represent  the  original 
function  after  it  has  passed  through  the  filters  corresponding  to  the  different  bands.  Figure  7  shows 
the  projections  for  a  sampled  chirp  function  and  analysing  Daubechies  wavelet  with  L=7.  As  is  to 
be  expected,  the  projections  migrate  across  the  frequency  plane  as  time  advances.  Figure  8  shows 
the  projections  for  a  cubed  sine  wave.  It  will  be  noted  that  the  signal  has  been  split  into  components 
that  exhibit  its  two  constituent  frequencies. 
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Consider  some  sine  wave  data  (Figure  9)  and  some  noisy  sine  wave  data  (Figure  10).  (Both  data 
sets  consist  of  128  samples).  Figures  11  and  12  show  the  respective  wavelet  expansion  coefficients 

for  the  Daubechies  wavelet  with  L=4.  Starting  at  the  bottom,  the  horizontal  rows  represent  J', 

t/  '.  t/  '.  t/ t/  '  and  d"*  (negative  values  indicated  by  the  lighter  shading).  Quantities  and 

d  '  are  represented  in  order  by  the  vertical  columns.  (The  number  in  the  top  right  hand  box  is  the 
figure  height  in  the  coefficient  scale).  As  is  to  be  expected,  the  noise  components  are  confmed  to 
the  finer  scales. 


6  SIGNAL  STRUCTURE 


The  wavelet  coefficients  { d‘ )  can  provide  valuable  information  concerning  the  regularity  of  a 
signal  [81,  regularity  described  in  terms  of  Lipschitz  exponents.  Function  /  is  Lipschitz  a  at  if 
there  exists  positive  A  and  such  that 

(36) 

for  some  polynomial  (order  n)  and 

Function  /  is  uniformly  Lipschitz  a  over  an  interval  if  (36)  holds  for  all  pairs  of  points  (x.^q)  in 
that  interval.  The  following  result  [9]  can  be  faroved: 

For  a  wavelet  with  greater  than  n  vanishing  moments,  a  function  /  is  uniformly  Lipschitz  a  if  and 
only  if  there  exists  K>0  such  that 


|d/|  <  K2^ 


(37) 


for  all  j  . 

If  the  wavelet  has  compact  support,  it  is  obvious  that  the  coefficients  at  finer  scales  will  provide 
localised  information  concerning  the  Lipschitz  exponents  a  of  the  function. 

Although  the  above  result  does  not  extend  directly  to  negative  exponents,  the  behaviour  of  the 
wavelet  coefficients  can  still  provide  a  good  indication  of  a.  Figure  13  shows  data  sampled  from  a 
function  that  contains  a  pulse,  a  step  and  a  continuous  transition.  The  corresponding  wavelet 
coefficients  (Daubechies  wavelet  with  L  =  10)  are  shown  in  Figure  14.  It  will  be  noted  that  the 

coefficient  maxima  exhibit  the  expected  behaviour  (the  features  have  exponents  a  with  values  -1,  0 
and  1  respectively).  This  ability  to  provide  information  about  transient  features  constitutes  one  of  the 
major  attractions  of  a  wavelet  approach. 
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Local  results  are  more  difficult  to  denve,  but  the  following  theorem  has  been  proved  by  Jaffard  [9], 
If  /  IS  Lipschitz  a  at  x,^. 


\d/\  <C2~^  (Ul2Lt„-tr) 


(38) 


Conversely,  if  equation  (38)  holds  and  f  is  uniformly  C®  for  a  positive  number  P,  there  exists  a 
polynomial  P  (depending  on  .r,,)  of  degree  less  than  a  such  that 


f 

\f(.x)  -  P(x-x,)\  <  C|x:-x,J“  log. 


2 


jr-x„ 


(39) 


and  this  result  is  optimal. 


7  CONCLUSIONS 

The  present  report  has  described  discrete  wavelet  analysis  through  an  approach  that  was  pioneered 
by  Daubechies  [1,6].  This  approach  has  been  developed  into  a  computer  software  package 
(Appendix  B)  that  can  be  used  to  analyse  sampled  data.  The  report  has  included  many  examples 
that  were  produced  using  this  package  and  these  amply  demonstrate  the  potendai  of  a  wavelet 
apfx-oach  to  signal  analysis.  A  major  advantage  of  the  approach  is  its  ability  to  provide  a  degree  of 
time  localisation  and  this,  together  with  the  existence  of  very  fast  algorithms,  makes  it  a  powerful 
tool  for  analysing  signals  with  transient  features.  In  particular,  wavelets  have  found  an  important 
application  in  the  analysis  of  speech  [7],  It  should  be  noted,  however,  that  the  approach  only 
provides  spectral  information  in  terms  of  frequency  bins  that  have  constant  Q  and  so  may  not  be 
suitable  for  applications  requiring  a  high  degree  of  frequency  resolution. 
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V  V 


Figure  7  Decomposition  of  a  chirp. 
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Figure  8  Decomposition  of  a  cubed  sine  wave. 
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APPENDIX  A  COMPUTER  SUBROUTINES 


The  ft>llowing  FORTRAN  subroutines  are  the  minimuin  that  is  required  to  implement  a  wavelet 
analysis.  Coefficients  Cq  to  c,,.,  (all  other  coefficients  zero)  should  be  placed  in  the  array  elements 
c(  1)  to  c(N).  Subroutine  scale  places  the  values  of  scale  function  ^  at  points  0  to  N-1  into  the  array 
elements  pbid)  to  phi(N).  The  function  of  interest  is  sampled  at  M  points  (M  =  2*^)  with 
corresponding  values  f(l)  to  f(M).  Subroutine  fit  will  estimate  the  filtered  samples  that  are  used  in 
the  initial  approximation  \{  and  place  them  in  the  array  elements  a(l)  to  a(M).  The  wavelet 
analysis  is  performed  by  subroudne  dccomp  which  replaces  a(l)  to  a(M)  with  the  wavelet 

coefficients  . d').  Finally,  subr^tine  neons  can  be  used  to  undo  the  above  analysis 

to  a  level  where  there  are  ML  scaling  function  basis  terms. 


***  Values  of  scaling  function  phi  at  points  0  to  N-1  *** 

subroutine  scale (c,N, phi/ 

dimension  c (N) ,phi (N) ,a(51, 51) ,e(51) 

nl=N-l 

do  11  i=2,nl 
do  11  j=2,nl 
cc=0 . 0 

ij=2* { i-1) -j+l 

if (ij .ge.O.and.ij .It.N)  cc=c(ij+l) 

11  a( i-1, j -1) =cc 
em=0 .0000001 
n2snl-l 
n3=:51 

call  niaxev(a.e.em,n2,n3) 
phi(l) =0.0 
phi(M) =0.0 
do  22  i=2,nl 
22  phi(i)=e(i-l) 

ph=0.0 
do  33  i=l,N 
3  3  ph=ph+phi  ( i ) 

do  44  i=l,N 

44  phi  ( i) =phi (i) /ph 

return 
end 


•**  eigenvector  corresponding  to  the  maximum  eigenvalue  *** 


subroutine  m2ucev(a.e,em2tx,n,m) 

dimension  a(m,n) ,e(n) ,v(200) 

iffn.eq.l)  then 

e(l)=l. 

emax=a (1,1) 

return 

endif 

error=einax 
do  1  i=l,n 

1  e(i)=1.0 

99  do  4  i=l,n 

v( i) =0 .0 
do  4  j=l,n 

4  v(i)=v(i)+a(i, j) *e{j) 

eoldsemax 
do  2  i=l,n 

2  if (abs (v( i) ) .gt .abs (emax) )  emaxsv(i) 
do  3  i=l,n 

3  e ( i) =v( i) /emax 

if (abs ( (emax-eold) /eold) .It. error)  return 
go  to  99 
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end 


**’  production  of  filtered  s^unples  from  function  samples  **• 


subroutine  f it (a,M, f ,phi,N) 

real*4  a(m} , f (ns) .phi(N) 

do  1  i=l,ns 

a( i>  =0 .0 

do  1  j=l.N 

tx=l-*- j-2 

ix=mod( ix,M) 

i  a; i ) =a( i ) +phi ( j > *f (ix+1) 
return 
end 


convolution  of  arrays  x  and  g  *•* 


real»4  function  conv(g,n,na,x,m, 1) 

real*4  g (n) , x(m) 

convsO .0 

do  1  isl,n 

ix=l-i-na 

if(ix.lt.O)  ix=ix-*-m 
ix=mod( ix,m) 

1  conv=conv+g(i) •xCix+l) 
return 
end 


•••  wavelet  coefficients  from  filtered  samples  *** 
subroutine  decomp ( a, H,c,N) 

real*4  a(M) ,c(n) ,gl(100) ,g2(100) , ft (1024) ,dt(1024) 

nl=-N+l 

n2=-l 

sqr2=sqrt (2 . ) 
do  5  i=l,N 
gl{i)=c(N-i+l) /sqr2 
5  g2(i)=c(i)*real( (-l)**i)/sqr2 

nt=M 

1  nn=nt 
nt=nt/2 

do  2  j=l,nt 
k=j*2-l 

2  ft( j)=conv(gl,N,nl,a,nn,k) 
do  3  j=l,nt 

k=j»2-l 

3  ft{ j)conv(g2,N,n2,a,nn,k) 
do  4  j=l,nt 
a(j+nt)=dt(j) 

4  a(j)=ft(j) 

if (nt.lt. 2)  then 

return 

else 

go  to  1 

endif 

end 


*•*  reconstitution  of  s2unples  from  wavelet  coefficients  *** 
subroutine  recons  (f,ML,c,N) 

real*4  f (M) ,c(n) ,gl(100) ,g2{100) , ft(1024) ,dt(1024) 

nl=0 

n2=-N+2 

sqr2ssqrt(2.) 

do  5  i=l,N 

gl(i)sr(i) /sqr2 
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5  g2(i)=real(  +  )  *c(N-i-*-l) /sqr2 

nt  =  l 

1  nn=nt 

do  2  j=l,nn 
ft(2*j-l)=f (j) 
ft(2*j)=0.0 
dt(2*j-l)=f (nn+j) 

2  dt(2*j)=0,0 
nt=nt*2 

do  3  j=l,nt 

3  f  ( j  )  =conv(gl.N,nl,  f  t ,  nt .  j )  •♦■conv(g2, N.  n2, dt , nt .  j ) 
if(nt.ge.ML)  then 

return 
else 
go  to  1 
end  if 
end 
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APPENDIX  B  SOFTWARE 


The  WV  wavelet  software  package  runs  on  a  DOS  PC  and  has  the  following  facilities: 

1)  To  choose  a  Daubechies  wavelet  and  to  display  the  wavelet  function,  the  scaling 
function  and  their  respective  Fourier  transforms. 

2)  To  enter  and  filter  data. 

3)  To  display  data. 

4)  To  display  a  frequency  analysis  of  data  (the  amplitude  of  the  discrete  Fourier 
transform). 

5)  To  change  to  a  logarithmic  data  scale. 

6)  To  change  to  an  exponential  data  scale  (reverses  5). 

7)  To  calculate  the  wavelet  coefficients  corresponding  to  user  supplied  data.  The 

coefficients  are  displayed  in  the  format  of  Figures  11,  12  and  14.  For  a  data  set  of 
initial  size  2*^  (data  sampled  at  unit  intervals),  the  rows  from  bottom  to  top 
represent  the  coefficients  of  projections  D_,/,D.j/, ..,  The  first  column 

represents  the  coefficients  of  A_j^/  and  tfie  second  column  represents  the 
coefficients  of  D.^/.  The  top  right  hand  box  shows  the  display  height  in  the 
coefficient  scale. 

8)  To  display  the  projections  (only  those  projections  above 

a  certain  threshold  will  be  displayed). 

9)  To  output  a  graphics  screen  as  a  file  of  postscript  instructions. 

10)  To  exit  the  package. 

It  should  be  noted  that  a  graphics  screen  is  exited  by  pressing  the  enter  key  and  that  the 

software  can  accept  at  most  1024  data  samples  (free  format). 

The  computer  disc  contains  the  WV  software  and  some  sample  data. 
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